Kelvin-Helmholtz instability in two-component Bose gases on a lattice 
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We explore the stability of the interface between two phase-separated Bose gases in relative 
motion on a lattice. Gross-Pitaevskii-Bogoliubov theory and the Gutzwiller ansatz are employed 
to study the short- and long-time stability properties. The underlying lattice introduces effects 
of discreteness, broken spatial symmetry, and strong correlations, all three of which are seen to 
have considerable qualitative effects on the Kelvin-Helmholtz instability. Discreteness is found to 
stabilize low flow velocities, because of the finite energy associated with displacing the interface. 
Broken spatial symmetry introduces a dependence not only on the relative flow velocity, but on the 
absolute velocities. Strong correlations close to a Mott transition will stop the Kelvin-Helmholtz 
instability from affecting the bulk density and creating turbulence; instead, the instability will excite 
vortices with Mott-insulator filled cores. 
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I. INTRODUCTION 

The Kelvin-Helmholtz instability is a dynamical insta- 
bility of the interface between two fluids that move rela- 
tive to one another. This hydrodynamical instability can 
occur in very different settings, for example when the 
wind is blowing over the surface of the smooth ocean 
or in plasma flows in the Earth's magnetic field 1|, 
or at the AB-phase boundary in superfluid "^He 0, 3|. 
The creation of degenerate mixtures of weakly interact- 
ing different bosonic gases enables studies of the 
Kelvin-Helmholtz instability in quantum systems which 
are more controllable as well as easier to measure in de- 
tail than more strongly interacting superfluids such as 
liquid ^He. Takeuchi, Suzuki et al. 0, [1] performed 
theoretical studies of the Kelvin-Helmholtz instability 
in a two-component condensate. Essentially, Ref. 
confirmed the expectation that Bose-Einstein condensed 
atomic gases provide an ideal setting to study the clas- 
sic Kelvin-Helmholtz instability without viscosity com- 
plicating the picture. Ref. [8] found a considerably al- 
tered instability dispersion relation in the case of a wide 
interface with large overlap between the two condensates; 
this was classified as a counter-superflow instability. 

Putting a Bose gas in an optical lattice introduces sev- 
eral new features. Most dramatically, the gas exhibits 
a quantum phase transition between a superfluid and a 
Mott insulating state when the ratio between the tun- 
neling and interaction energies passes a critical value Q . 
In a binary condensate the phase diagram is richer, dis- 
playing different combinations of Mott insulator and su- 
perfluid in coexisting or phase separated configurations 
[Tol[Tl|. Moreover, even in the superfluid state where a 
significant portion of the atoms are Bose-Einstein con- 
densed, the discreteness and the broken translational 
symmetry can have decisive effect on the motion of the 
bosons, as we shall see. 

In this paper, we study the Kelvin-Helmholtz-type in- 
stabilities of a phase separated binary condensate on a 
lattice. In Sec. [TTl we present the equations that we work 



with. In Sec. IIII Al we study a weakly segregated pair 
of condensates in order to explore the effects of broken 
translational symmetry on the interface dynamics. In 
Sec. IIIIBl we study instead a strongly segregated pair of 
condensates to see the effects of discreteness of the lat- 
tice. Sec. IIVI presents a heuristic derivation of a closed- 
form expression for the numerically obtained instabilities. 
In Sec. |Vl we study the system close to a Mott transi- 
tion in order to explore the effects of strong correlations. 
Finally, in Sec. IVIi we summarize and conclude. 



II. THEORY OF A TWO-COMPONENT 
LATTICE BOSE GAS 

The system is assumed to consist of two species of 
bosons at nearly zero temperature hopping on a square 
lattice, as can be realized with evaporatively cooled 
atoms in an optical lattice [4-6]. If the lattice is deep 
enough, the tight-binding approximation is valid and the 
bosons can be assumed to occupy the lowest band only. 
For simplicity, we assume the two species to have equal 
masses, interaction constant and tunneling properties; 
this can be realized by choosing the two components to 
be two spin states of the same element. The many-body 
Hamiltonian governing this system is 
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Here, Jj are the tunneling matrix elements 
in-species interaction parameters and U12 is the inter 
species interaction parameter, is the chemical poten- 
tial of species j . The summation index j = 1,2 denotes 
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the different species and tlie index r runs over tfie lattice 
sites. We will be studying a two-dimensional square lat- 
tice in this paper, and we therefore pass to a vector nota- 
tion, r = (x, 2:), where we let the dimensionless Cartesian 
coordinates x and z take on integer values. 

For strong enough hopping and for high enough den- 
sity, the lattice gases are almost entirely Bose-Einstein 
condensed and the problem can be treated using Gross- 



Pitaevskii and Bogoliubov analysis. This regime is suit- 
able for identifying effects of the discreteness and broken 
translational invariance of the lattice, effects that do not 
depend on quantum fluctuations. In this regime, the lat- 
tice gas is accurately described using condensate wave- 
functions $j(r,t) — {cijr), whose dynamics follows from 
the discrete two-component Gross-Pitaevskii (GP) equa- 
tion 



= -JiV2<J>i(r,t) + [C/n|<J>i|2 + C/i2|<I>2|2] $i(r,t). 



" dt =-'^2V2$2(r,0+[C^22|$2|' + C/i2|$iP]$2(r,t), (2) 

I 

where the discrete Laplacian is defined as dependent perturbation, 

V^$(r) = ^$(r'), (3) 

$i(r,t) = [iii{z,t) + S^i{x, z,t)]exp{-ifit), 

and the sum over r' runs over the nearest neighbors to ^ / J-^ r,Tf r +\ 1 ;:it, r j^m / • .^^ fA\ 

^j^g gj^g J. ^ *2(r, t) = [^'2(2;, t) + 5*2(2;, z, t)] exp{-ifit). (4) 



A. Bogoliubov approach 

The equation of motion to zeroth order 
In the Bogoliubov approximation one assumes for each independent ID GP equation that is used 
component a stationary wavefunction with a small time- the interface profile, 
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is the time- 
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[C/ll|*l(z)|2 + C/l2|*2(z)P] *l(^) - Ml'J'l(z), 
[C/22|*2(2)P + C/l2|*l(z)P] ^2{Z) = /^2*2(2), 



(5) 



where we introduced notation analogous to Eq. (jS]), 
92 ^-^.(z) 



where 
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^>,{z + l) + ^,{z-l) 



(6) 



Anticipating that the nonlinear term in the GP equa- 
tion will couple positive- and negative-frequency modes, 
we define in a standard way 

(5*i(r,t) = [ui(r, t) exp(— iwt) + w,*(r, exp(zajt)] exp(— i^i). 
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Ignoring the second order terms in the perturbations ^^'^ 
we get the Bogoliubov equations for the discrete double 
condensate system. 



Mv — ujfjv, 



(7) 



U2{r) 



(10) 
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The coefficients are defined tlirougli C = ?7i2^'i^'2! ^ 
C/i2*i*2, and Ui = l^-ip, 



i7i = -JiV2 + 2C/iini + i7i2n2-Mi, (H) 



and 



i?2 = - J2V2 + 2C/22n2 + J7i2ni - H2. (12) 



witli eigenvalue w. Combining Eqs. (jl4m5|) to obtain 



Let us repeat some common terminology and standard 
results [l^: The norm of a Bogoliubov mode is defined 
as 



^v\v)f]w{Y) 



(13) 



Unless the norm is zero, we are free to rescale the Bo- 
goliubov vector so that its absolute value is unity. A 
mode with ||v|| = 1 is called a quasiparticle mode and 
a mode with ||v|| = —1 is called a quasihole mode. If 
^ is the ground state, then all the quasiparticle modes 
have positive eigenfrequencies and all quasihole modes 
have negative eigenfrequencies; otherwise, no such rule 
applies. A mode with complex eigenfrequency has zero 
norm. Such modes will grow exponentially in time and 
thus they signal a dynamical instability. We will there- 
fore be interested in the imaginary parts of the Bogoli- 
ubov eigenfrequencies in this article. 

The symmetries of the problem give rise to degen- 
eracies in the Bogoliubov spectrum. First, the well- 
known quasiparticle-quasihole symmetry of the Bogoli- 
ubov equations [l2| carries over to the two-component 
case: Each real eigenvalue w with an eigenvector v = 
(Mi(r), i>i(r), U2(r), W2(r))"^ has a corresponding solution 
with an eigenvalue — w and an eigenvector 



V = - (vi{v)\m{v)\v2{v)\u2{vrY . (14) 

In addition, imaginary eigenvalues appear in pairs of lj 
and uj* . Since we will be mostly interested in a symmetric 
interface where ^'i(r) — 5'2(— r), we have an additional 
degenerate solution 



{V*2{-V),ul{-V)y,{-V),ul{-V)Y 



(16) 



with eigenvalue —uj. The latter two symmetries hold 
strictly only when equal and opposite currents are in- 
duced in the two components, as we will discuss further 
below. 



1. One moving condensate 

The discrete system is not translationally invariant, 
and we can therefore not use the relative velocity of the 
condensates as a unique parameter; also the absolute in- 
dividual currents are physically relevant. In particular, 
we will study two different cases: the case where only one 
of the condensates carries a current, and that where the 
condensates support equal but counter-flowing currents. 

We first assume that component 2 is moving along the 
X direction and that all other position dependence of the 
stationary solutions is along z. Furthermore, the sta- 
tionary wavefunctions can be taken as real. This means 
'^i{y) = *i(z) and \I'2(r) = ^2(^:) exp(ifca;), with -tt < 
k < TT. We then have C = Ui2^i{z)'^2{z)e^p{-ikx), 
D = [/i2*i(z)*2(2)exp(ifca;), and 



-J2V^+2U22n2 + Ui2ni-H2-2J2{l-cosk), (17) 



where /i is the chemical potential of the second compo- 
nent in the absence of a current; the term 2J(1 — cosfc) 
due to the current is made explicit for convenience. The 
excitations can be expanded as plane waves along the x 
direction. Dependency on the x coordinate is eliminated 
with the choice 

Mi(r) ~ ui{z) exp{ikeX — ikx/2), 
vi{r) = vi{z) exp{ikeX — ikx/2), 
M2(r) — U2{z) exp{ikeX + ikx/2), 
V2{r) — V2{z) exp{ikeX — 3ikx/2). (18) 

The Bogoliubov equations now reduce to a ID problem, 
and the Bogoliubov matrix takes the form 



M- 



ID 




(19) 



(u2(-r),W2(-r),'Ui(-r),z;i(-r)) , 



(15) where the diagonal elements are 



Hiu = Hiv = - JiTTT + [1 - cos{ke - k/2)] + 2Uiini + Ui2n2 - Mi, 
oz^ 

H2.U2 = -J2^ - 2J2(1 - COsfc) + 2J2 [1 - COs{ke + k/2)] + 2U22n2 + Ui2ni - 



2,1)2 



■dz^ 



-2J2(1 - cos A:) + 2J2 [1 - cos(/ce - 3fc/2))] + 2U22n2 + Ui2ni - /X2. 



(20) 
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2.. Counter- directed currents 



In the case of equal and opposite currents in the re- 
spective components, we write ^"1(0;, z) = ^'i(z)e*'^^/^ 
and '^2{x,z) = '^2{z)e~^^^/'^ . The excitations can now 
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Expressed in the reduced ID Bogohubov amphtudes 
Uj{z),Vj{z), the symmetries (|14m6p acquire new forms, 
since the effective Hamiltonians appearing in the diagonal 
of the Bogoliubov eigenvalue problem are different. The 
symmetry of greatest interest to us is that for a given k 
and ke , each real eigenvalue w with a reduced eigenvector 
^1(2:), U2{z), V2{z))'^ now has a corresponding hole 
mode with eigenvalue —uj and reduced eigenvector 

{v2{-z),U2{~z),vi{-z),ui{-z))'^ . (23) 

III. KELVIN-HELMHOLTZ INSTABILITY IN A 
DISCRETE SYSTEM 

After assuming that the two condensates have equal 
chemical potentials /zi = /i2 = tunneling Ji = J2 = J, 
and interaction strengths Un = U22 = U, the system can 
be characterized by three parameters. These are chosen 
to be the ratios J/U and U12/U, and the number density 
per site far from the interface, n. The chemical potential 
fj, can then be calculated for given J, U12, and n. In ad- 
dition, the problem is characterized by the wave vectors 
for the flow in each of the components, (fci, ^2)- In Sec. 
im we derived the detailed form of the Bogoliubov equa- 
tions for one moving condensate, {ki,k2) = (0, — fc); and 
opposite currents, (fci,fc2) = {k/2,—k/2). 

A. Weak coupling and broken translational 
symmetry 

We first consider the case J/U — 1, U12/U ~ 1.1, 
n = 1. Figure [T] shows the stability properties of the 
system in the case of one moving condensate, (fci, k2) = 



be written 

ui(r) =ui(z)e''==^+''=^/2^ 
wi(r) =^;l(z)e*'=="-*'="/^ 

V2{v) =V2{z)e'^'''+'^''''^. (21) 

We find that the off-diagonal terms of the matrix Mid 
are the same as before, but to the diagonal we get 



- Ml, 

- Ail, 

- M2, 

(22) 
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FIG. 1. (Color online) Growth rate of the dynamically most 
unstable mode (i.e. the imaginary part of the mode with the 
largest imaginary frequency) as a function of flow wavevector 
k and and mode wavevector ke (note that a shift fee — fc/2 has 
been applied to the y-axis to make the symmetry around fee = 
k/2 more transparent). Here, one component is in motion 
with wavevector k and the other is stationary. The lattice has 
256 sites, the parameter values are J/U = 1, Uu/U = 1.1, 
and the asymptotic density is n = 1. 



(0, — fc). The figure plots the largest imaginary excitation 
frequency as a function of flow wavevector fc and excita- 
tion wavevector kg,; thus, a vertical cross-section in the 
figure at a fixed fc gives the imaginary part of the spec- 
trum for a given current. The lattice is of finite size, 256 
sites, in the numerics (in fact, real-life lattices may be 
of similar dimensions) and this infiuences the numerical 



- cosfc/2) + 2Ji [1 

- cosfc/2) -f 2Ji [1 

- cosfc/2) -f 2J2 [1 

- cosfc/2) -f 2J2 [1 



- cos(fce + fc/2)] 

— COs(fce — fc/2)] 

— COs(fce — fc/2)] 

- COs(fce -I- fc/2)] 



+ 2[/iini -I- Ui2n2 
+ 2C/iini -I- 7/12^2 
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+ 2U22n2 + £^12^1 



5 



value of the imaginary part somewhat. The effect is such 
that the imaginary part becomes larger as the lattice size 
is reduced to 128. The physics at low wavevector k is 
clearly recognized from the case of a pair of classical flu- 
ids or a homogeneous two-component BEC ^, i8|, lH, [l3| • 
For a given fc, not too large, a range of excitations with 
wavevectors < \ke — k/2\ < k become unstable. The 
amplitudes of the unstable modes are localized at the in- 
terface, indicating that these are indeed interface modes. 

We compare to the case of a continuous, i.e., non- 
lattice system. In Ref. 8], it was observed that in the 
case of a narrow interface, the binary condensate displays 
the classical dispersion relation [13| 

u = ^J^^k^~vVt^ (classically), (24) 

where v is the velocity of the unperturbed flow and a is 
the surface tension. However, if the interface is wider, the 
upper stability line lies close to the line \ke — k/2\ = k; 
this was in Ref. Q termed the counter-superflow insta- 
bility. Apparently, Fig. [1] is in this latter regime, and the 
discreteness is not seen to cause considerable deviations 
from the linear stability line. 

We note a few curious features in Fig. [T] For very 
small k, the instability seems to be absent; this will be 
discussed further in Sec. 1111 Bl For wavevectors exceed- 
ing a critical value, A: > fee = 7r/2, it is seen in Fig. [T] 
that all excitation wavevectors become unstable. This is 
no longer an interface mode, but an instability known to 
occur also in single-component discrete BECs [l5|. The 
physics behind the instability is that the single-particle 
dispersion relation, = J(l — cosfc), has an inflection 
point at kc, and thus the effective mass goes from pos- 
itive, to diverging, to negative. This is well understood 
and studied in single-component BECs and we will not 
dwell on it further, but only conclude that the interesting 
range of k values for the case (^1,^2) = (0, — fc) ranges 
between — 7r/2 and tt/2. 

To illustrate the long time behavior of Kelvin- 
Helmholtz instability, we simulate the time development 
of the discrete GP equation. A series of snapshots are 
shown in Fig. [5] for the case k = 57r/64 « 0.245. The 
time development is simulated using a split-step Fourier 
technique, with a grid of size 64x64 sites and periodic 
boundary conditions along both directions. For the cho- 
sen wavevector k — 57r/64, a surface wave with a domi- 
nant wavevector fee = 37r/64 can be seen to grow expo- 
nentially until the nonlinear stage is reached, where the 
condensates break up into smaller fragments. 

For the case of counterflow, (fc/2, — fc/2), the result for 
the imaginary parts of the modes is shown in Fig. |3l By 
definition of k, the critical relative wavevector for neg- 
ative effective mass is in this case kc = tt. For large 
enough A;, it is seen that the modes with small wavevec- 
tors ke are actually stabilized. We thus define a second 
critical wavenumber for the flow, ks, above which long- 
wavelength excitations are stabilized. The corresponding 
flow velocity is denoted . From Figure [31 we read off 




FIG. 2. (Color online) Time development of the instability 
for the case k = 57r/64, and remaining parameter values as 
in Fig. [1] Panels (a),(c),(e),(g) show the density of the 1 
component, where red indicates high density and blue is zero 
density; panels (b), (d), (f), (h) show the corresponding phase. 
Snapshots are taken at times t — 44 for (a)-(b); t = 176 for 
(c)-(d); t = 265 for (e)-(f); and t = 400 for (g)-(h). 




0.5 1 1.5 2 

k 



FIG. 3. (Color online) Largest imaginary part of the exci- 
tation frequency as in Fig. [T] but with counter-flowing con- 
densates having wavevectors (fc/2, —fc/2). The lattice has 256 
sites, parameter values are J/U = 1, U12/U = 1.1, and the 
asymptotic density is fi = 1. 
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FIG. 4. (Color online) Comparison of the largest imaginary 
part of the excitation frequency as a function of the excitation 
wavenumber ke, in the cases of one flowing component, {k,0) 
(dashed red), and counterflow, (fc/2, — fc/2) (solid blue). Pa- 
rameter values are as in Fig. [Jl but the wavevector is fixed at 
= 1.3 in (a), and = 1.0 in (b). For the case of one flowing 
component, the ke axis is shifted by ke — ke — k/2. The lattice 
has 256 sites, parameter values are J/U = 1, U12/U — 1.1, 
and the asymptotic density is n = 1. 

kg ~ 1.5 (the closeness to 7r/2 is suggestive but, in fact, 
a coincidence). This is a counterpart to the classically 
known stabilization of long wavelengths when the flow 
velocity exceeds -\/8 times the sound speed c , 

Vs = VSc (classically). (25) 

However, in the discrete case, there seems to be no sim- 
ple relationship between the onset of long-wavelength 
stability and sound speed. In this particular case, the 
flow velocity in each component is Vs — 2Jsin(fcs/2) sa 
lAU in dim ensionless units, whereas the sound speed is 
c = V JUn = U, but further numerical experimentation 
shows only a weak dependence of Vs on c; this puzzle will 
be resolved in Sec. IIVI 

Since the system lacks Galilean invariance, the spec- 
trum for a given flow wavevector k is not invariant under 
a transformation (fci, fci) — )■ {ki — q, k2 + q). We plot as an 
example in Fig. 2] the imaginary frequency for the same 
parameters as above, for two choices of k, comparing the 
cases (fc, 0) and (fc/2, —k/2). It is seen that the curves are 
similar, but the symmetric case {k/2, —k/2) is more un- 
stable, in the sense that it has a higher maximum imagi- 
nary frequency. Further numerical experimentation (not 
shown here) indicates that this observation holds quite 
generally. 

B. Narrow interface physics: Effects of discreteness 

Figure [5] shows the stability spectrum of a system at a 
higher density, n = 4, and a stronger repulsion between 
the components, C/12 = 1.5C/. The combined effect of 
increased density, which means larger interaction energy. 



and increased inter-species repulsion, implies a narrower 
interface Already in the case n = l,C/i2 = 1.1J7, 

we saw that low flow velocities appeared to be stable, 
and here the effect is more pronounced. This feature of 
the instability is not seen in simulations of homogeneous 
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k 

FIG. 5. (Color online) Growth rate of the dynamically most 
unstable mode (i.e. the imaginary part of the mode with the 
largest imaginary frequency) as a function of flow wavevector 
k and and mode wavevector ke . Here, both components are in 
motion with counter-directed wavevectors {k/2, —k/2). The 
parameters are J/U = 1, U12/U — 1.5, and the asymptotic 
density is n = 4. 

BECs Hi Bi D or in classical fluids, in the absence of exter- 
nal forces We attribute this to the discrete density 
profile at the interface, as will be detailed below. 

In the figure, we note also that there is a stable region 
for large flow wavevectors k > ks, just as in the case 
n = 1, but now ks « 2.8. We will return to this in Sec. 
ITVl Finally, we observe that the upper instability line at 
large k values is still approximately linear with unit slope, 
indicating that we are in the counter-superflow regime 
despite the narrower interface. 

We now turn to the stability of slow flow wavevectors 
k. First recall how a dynamical instability in the Bogoli- 
ubov equations result from a collision between a quasi- 
particle mode and a quasihole mode [13] • Starting from 
a reference point in phase space - here, take A; = - we 
consider a quasiparticle mode and a quasihole mode 
Vb with real, nearly degenerate eigenfrequencies ujQa and 
■^oh- We then slightly increase a control parameter - in 
this case, k - and thus transform the Bogoliubov matrix 
M M' = M + SM. Writing v — aVa + bvb and solving 
the Bogoliubov equations for the new eigenfrequencies w, 
we find in the general case 



(26) 
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FIG. 6. Few Bogoliubov eigenvalues for a fixed excitation 
wavevector fee = 0.3, as functions of the flow wavevector fe. 
Parameters are as in Fig. [S] Upper panel: Real parts, lower 
panel: Imaginary parts. 



where StUa ~ v^SMva, Suii, = — v^(5A/vb, and X = 
vl^SMvi,. In other words, when the two modes become 
degenerate, the cross-term makes them unstable. 

In the present case, we consider two modes related by 
the symmetry of Eq. (|23p : thus, we have woa + <5wa — 
—LUQb — Sujb, and 



± 



(27) 



In a homogeneous, i.e., non-discrete system, the lowest- 
lying surface mode corresponds to a uniform translation 
of the boundary and is thus a Goldstonc mode; at fc = 
fee = it has zero energy, i.e., woa = in the limit ke — 0. 
With increasing ke, its energy increases from zero. A 
finite k results in a shift in the Bogoliubov matrix SM 
as discussed above and the Kelvin-Helmholtz instability 
comes about because |X| is always greater than Scda hi 
these systems. The onset of instability therefore tends 
towards k — as k^. tends to zero. 

In a discrete system, however, translational symmetry 
is lost and there is no Goldstone mode: A small shift 
of the boundary profile costs a finite amount of energy. 
This translates into woa = e > at /ce = 0, which in turn 
creates a threshold value of k for the onset of instability. 
For stronger interactions or higher densities, the inter- 
face is thinner, and the effects of discreteness are more 
pronounced. This explains the increased threshold value 
for higher density. 

Figure [S] displays some of the lowest-lying Bogoliubov 
eigenvalues for a fixed k^ = 0.3, as functions of fc, at n = 
4. It is clearly seen how instabilities result from a merging 
of two eigenstates with equal but opposite energies, thus 
related by the symmetry in Eq. (j23p . For large k similar 
instabilities happen for many other modes and here, for 
clarity, we plot just the appearance of the first instability. 
The modes responsible for the dominant instability start 
out at a small but finite energy for fc = 0; this energy 
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FIG. 7. Bogoliubov eigenstates for the lowest-lying sur- 
face modes, functions Uj{z) (blue solid lines) and Vj{z) (red 
dashed) at fee ~ 0.3 and fe = 1.25, which is below the instabil- 
ity threshold. Upper panels: Quasiparticle mode; Lower pan- 
els: Quasihole mode. Left panels: mode functions ui(2:), 111(2:) 
associated with component 1, and right panels: mode func- 
tions 112(2), V2{z). The lattice has 256 sites, parameter values 
are J/U = 1, U12/U = 1.5, and the asymptotic density is 
n = 4. 



does not tend to zero as fee —> (although the latter fact 
cannot be inferred from the figure) . The mode functions 
Ui,Vi are plotted in Fig. [71 These are the lowest- lying 
surface modes corresponding to translation of the surface; 
notice how the symmetry (j23p is manifest in the plots. 



IV. TOWARDS AN ANALYTICAL 
DESCRIPTION 

In order to obtain a feeling for the mathematics and 
physics behind the Kelvin-Helmholtz instability, we show 
here a heuristic derivation of an approximate analytical 
formula for the excitation frequencies. In effect, we re- 
place the spatially dependent coefficients in Eq. by c- 
numbers. This can be justified by expanding the solution 
V in a suitable basis 18]. It can also be effected by sim- 
ply assuming a local-density approximation at the center 
of the interface. We discuss the more symmetric case 
with counter-directed currents, (fci,fc2) = (fc/2,— fc/2), 
since it allows for solutions on closed form. The matrix 
is approximated as 



M 



/e + A+ 

9 
C 

V c 
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9 
C 

c 



c 

c 
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C 
C 
9 
A+ 
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where g is proportional to the in-species interaction en- 
ergy, and C is a measure of the inter-species repulsion, e 
is a self-energy associated with bending in the z direction 
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of the respective mode function, which is to be identified 
with the finite energy ujQa of the almost- Goldstone modes 
dicussed above. Finally, A± are the relative kinetic en- 
ergies. In the discrete lattice case. 



A± = -2 J(l - cos A:/2) + 2 J [1 - cos(fce ± k/2)] , (28) 



This equation is of the same form as that in Ref. [8|, 
describing the counter-superflow instability. Clearly, the 
case with a minus sign lets become negative for certain 
parameter values, and thus describes the unstable modes. 
The discrete case yields somewhat clumsier expressions, 
but can be solved; we discuss the differences below. 

The choice of parameters can be done with the benefit 
of hindsight; however, assuming a local-density approxi- 
mation in the midpoint of the interface, where in the limit 
of a broad interface the densities ni and n2 are half their 
bulk values, we readily obtain g = Un/2 = /i/2, and 
C — Ui2n/2 = {Ui2/U)g. The self-energy e is smaller 
than the other energies, and so is C — g, so we first put 
e = and C = g in Eq. ((30)) and find the points where 
turns negative; these are the boundaries for instability. 
We obtain the solutions 

fee = 0, 

fc2 - fc2 - 4/^ = 0. (31) 
Linearizing around the first instability, we obtain 

^' = lkl(e~lk'^+0{kt), (32) 

which is the well-known linearly increasing imaginary fre- 
quency, existing for fc^ < 4/i. The modes are stabilized 
again when k,, > k. The upper stability limit lies at 
ke = k as observed in Ref. However, when fc^ > 4^, 
small wavevectors fee are stabilized. This is similar to 
the well known classical stabilization at relative flow ve- 
locities exceeding -v/S times the sound velocity, as we dis- 
cussed in Sec. IIII Al 14] . However, even in the continuous 
case, the prefactor does not match with the classical re- 
sult (nor with Ref. 0); we simply attribute this to the 
quite different kind of dispersion relation. 

A small but finite e will introduce small corrections 
to the above. In order to salvage the linear dispersion at 
small excitation wavevectors ke — which is observed to be 
the case numerically - we must put C = g + e/2. With 
a finite e, slow fiow velocities are stabilized. This hap- 
pens only in the discrete case, but in order to bring out 



and in the continuous case without a lattice, 

A± = —{kl±kek). (29) 
2m 

The continuous case is simplest: Inserting the dispersion, 
putting h = m = 1 for convenience, we obtain the eigen- 
frequencies 



(30) 

I 

the qualitative features while displaying reasonably com- 
pact formulas, we show here the bastardized expression 
resulting from using a continuous-case dispersion relation 
together with a finite self-energy. Expanding the eigen- 
frequency in powers of k and fcg, one obtains 

- lkUlk',+e)-^k^kl + 0{k\kt,k^kt,kX), (33) 

which is positive for k^ < + 2e. Thus, slow flow ve- 
locities are seen to be absolutely stable. Classically, this 
only happens in the presence of gravity or an equivalent 
asymmetric force 0, [l3| , but in the discrete system, there 
is stabilization even in the absence of forces. 

The discrete-case dispersion relation leads to more 
complicated algebra, but all the qualitative conclusions 
remain the same. When e = and g = C, the frequency 
is again zero for k^ — and k^ — k; a finite e stabilizes 
small k. The most important modification is that of the 
sound-speed limit, where we now obtain stability at small 
kf. a k > kg, with kg given by 

2Jcos(^) = -^l/2 + v/mV4 + 4J2, (34) 

whose solution is always real. The continuous-case ex- 
pression is recovered when /i/J is small; then the flow 
velocity is = 2 J sinks ~ \/4/^: similar to the continu- 
ous case. 

Furthermore, for k > 7r/2, the spectrum changes char- 
acter completely and a wide range of excitation wavevec- 
tors fee become unstable; this is again a manifestation 
of the bulk instability associated with negative effective 
mass of Ref. [lH| ■ 

Thus, a heuristic constant-amplitude calculation has 
managed to bring out all the features present in the nu- 
merical results. The familiar linear increase of imaginary 
frequency at small fee for a given flow wavevector fc; the 
stabilization at excitation wavevectors fee exceeding fc; 
the absolute stabilization of slow relative velocities; the 
stabilization of small k^ for k above a cutoff depending on 
the sound velocity; and even the bulk instability when the 
effective mass becomes negative, could all be described. 



I 

= {2g + hi + e){e + hi) + hlk' ± 2^e+ h^^j {2g + . + hl)hlk^ + C^{e + h^). 



9 



1 

0.9 
0.8 
0.7 


Psep Ml ) 






Psep SF 


MI+SF it 






0.6 








0.5 








0.4 






Coex SF 


0.3 








0.2 


Coex Ml f 






0.1 








n 









0.02 0.04 0.06 0.08 0.1 

J/U 



FIG. 8. Phase diagram for a mixture of two identical Bose 
gases, calculated in the Gutzwiller approximation. The chem- 
ical potential is taken to be /.i = 0.7U. The points A and B 
mark the parameters for which the two simulations described 
in the text are done. 

V. STRONGLY CORRELATED REGIME 

The discussion so far has been hniited to the case of two 
Bose gases in the Bose-Einstein condensed state. This 
is the state of the system when the ratio J/U is not 
very small and when the density is not very small. In 
the general case, one must account for quantum fluctu- 
ations whose most striking effect is to drive a quantum 
phase transition between the superfluid and Mott insu- 
lating phases [9]. In addition, quantum fluctuations are 
seen to shift the transition line between mixed and phase- 
segregated phases, as we shall see below. 

The many-body Hamiltonian ^ has been seen to dis- 
play a rich phase diagram [13, [HI- In Figure HI we display 
the portion of greatest interest for our purposes. The 
phase diagram is calculated in the Gutzwiller approxi- 
mation, which is explained below. The general picture 
is that for weak enough hopping J, except in degener- 
ate cases, both gases are in the Mott insulating state 
(marked MI in Fig. (|8])) with suppressed on-site number 
fluctuations and zero superfluid density; conversely, for 
large enough J the bosons are always superfluid (SF). 
If the inter-component interaction U12 is large enough 
compared with the chemical potential /i, then the two 
components are phase separated (marked Psep), either 
Mott insulating or superfluid; and if U12 is small enough 
or fi big enough, the two components coexist (marked 
Coex). There exists one further phase close to the line 
U12 — fi at weak hopping, where one component is Mott 
insulating and the other one is superfluid with a lower 
density (marked SF-I-MI in the phase diagram). Finally, 
note that for a pair of Bose-Einstein condensed gases, the 
transition between phase separated and coexisting phases 
takes place at U12 = U. In contrast, in the limit of weak 
J, phase separation takes place at U12 = fi- 



The Gutzwiller approximation is based on a mean-field 
ansatz for the many-body state, 

2 

iV'G(o)-nni'^'-.^-w)- (35) 

r j=l 

This ansatz amounts to treating both the hopping be- 
tween sites and the inter-species interaction in a mean- 
field fashion. This can be done since on-site entanglement 
between the two species is not expected to be important. 
For computations, it is convenient to expand the on-site 
states in a local Fock basis with an upper cutoff rimax, 

\(brAt)) = E a,j,„(t)|n),.,,. (36) 

In the Gutzwiller approximation, the local total density 
n and condensate wave function z for each of the compo- 
nents is readily computed. We characterize the system 
by studying the behavior of local density Uj, the local 
condensate density Ucj = I'&jP, and phase ipj = arg<i>j. 
In order to calculate the ground state, the Hamiltonian 
([T]) is minimized with respect to the complex coefficients 
Cr.j,n- To simulate dynamics, we propagate the coupled 
equations of motion [l^) HO] 

'^-^-J^^^^o{t)\H\^a{t)). (37) 

For the present simulation of Kelvin-Helmholtz physics, 
the initial condition was constructed by calculating the 
phase separated ground state, then imprinting counter- 
directed currents on the two components, and finally 
adding a small amount of random noise onto each of the 
complex coefficients Crj",n- 

We simulate interface dynamics in the "Psep SF" 
phase. In the Mott insulating phase, dynamics is quite 
trivially absent. In the "SF4-MI" phase, there is also 
no interface dynamics: Since the Mott-superfluid phase 
transition is second order, there is no surface tension and 
hence no capillary waves are expected to form between 
superfluid and Mott insulating regions. Phenomena con- 
nected to melting of Mott insulators are conceivable at 
high enough energies (cf. [13,[ll|), but we will not pursue 
this here. 

The first simulation is done at the point marked A in 
the phase diagram of Fig H J = 0.045C/, U12 = 0.85C/, 
and [I — 0.7U. Here, we have two phase separated super- 
fluids, and as seen in Fig. HI the Bose-Einstein condensed 
part of the fluid is less than half the total density. Also 
the superfluid density is expected to be small close to 
the Mott transition 22]. Figure [TU] shows the dynamics 
of a system with symmetrically imprinted wavevectors 
(fc/2,-fc/2), where k = 7n/16 « 1.3744. It is seen that 
the condensate densities behave very much like in the 
pure-condensate case. The time scale is vastly different 
due to the difference in tunneling matrix element J. In 
addition, it is seen that the total density is affected very 
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FIG. 9. (Color online) Cross sections of the density pro- 
file of two phase segregated Bose gases at points A (upper 
panel) and B (lower panel) in the phase diagram, respectively. 
Squares denote the density for component 1 ni; crosses de- 
note Uci; asterisks 712, and circles nc2- Lines are to guide the 
eye. 




FIG. 10. (Color online) Time development in the Gutzwiller 
approximation of a system with J — 0.04517, U12 ~ 0.85(7, 
and fj, — 0.7U, with wave vectors (fc/2, —k/2) where k — 
77r/16. Snapshots displayed in the different rows are taken 
at times (from top to bottom) 100/i7, 200/C/, 300/?7, and 
4000/f/. Left column indicates total density n, middle col- 
umn condensate density ric, and right column the phase (f, 
of bosons of component 1. The corresponding quantities for 
component 2 are not shown, but are similar to the ones shown 
mirrored in the x axis. 



little, and on the time scales we have simulated the sys- 
tem seems to be settling down into a phase separated 
steady state. The two fluids do not mix; the main effect 
of the KH instability is to excite vortices within the re- 
spective components. In the phase plots one can see vor- 
tices as phase singularities; they show up as zeros of the 
condensate density but the total density is not depleted. 



FIG. 11. (Color online) Time development in the Gutzwiller 
approximation of a system with J — 0.03517, U12 = 0.85C/, 
and jj, — 0.7(7, with wave vectors {k/2, — fc/2) where k = 
77r/16. Snapshots displayed in the different rows are taken 
at times (from top to bottom) 100/17, 200/C/, 300/17, and 
4000/(7. Left column indicates total density n, middle col- 
umn condensate density ric, and right column the phase ip, 
of bosons of component 1. The corresponding quantities for 
component 2 are not shown, but are similar to the ones shown 
mirrored in the x axis. 



These are vortices filled with Mott insulating atoms, as 
previously described in Refs. (23ll23|. 

We next show a simulation at the point marked B in 
the phase diagram of Fig [i J = 0.035(7, U12 = 0.85(7, 
and ^ = 0.7(7. This is still in the regime of two phase 
separated superfluids, but even closer to the Mott tran- 
sition. Figure [TT] shows the dynamics. We see that the 
dynamics is in fact more rapid here than in the case of 
larger J. At times t = 200/(7 and t = 300/(7, a signifi- 
cant wave pattern develops in the condensate component, 
and a.t t — 1000/(7, the system contains many more vor- 
tices than in the case in Fig. [TOl The total density is less 
affected, but exhibits some protruding features at the 
surface. Clearly, with a smaller condensate fraction, the 
threshold energy is smaller for modulation of condensate 
density and formation of vortices. Close to the transition, 
where the condensate density decreases rapidly towards 
zero, this overcomes the effect of the modest reduction of 
J. 



VI. CONCLUSIONS 

In summary, we have studied the Kelvin-Helmholtz in- 
stability on the interface between two lattice Bose gases 
in relative motion. The instability is seen to be affected 
by three effects introduced by the lattice potential: bro- 
ken translational symmetry, discreteness, and quantum 
fluctuations. Broken translational symmetry affects the 
instability in such a way that the excitation frequencies 
do not depend only on the relative velocity, but on the 
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flow velocities of the two gases separately; the symmetric 
case of two counter-flowing gases is seen to be the most 
unstable situation. Second, the discreteness will stabilize 
slow relative currents, so that the instability is prevented 
if the relative velocity is low enough. Finally, strongly 
correlated physics affects the physics in such a way that 
a neighboring Mott insulating phase will prevent the two 
Bose systems from mixing after the Kelvin-Helmholtz in- 
stability is first excited. Close to the Mott phase tran- 
sition, only the superfluid density will take part in the 



instability, but the total density will hardly be affected. 
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